library(brms)
library(data.table)
library(latex2exp)

# SET WORKING DIRECTORY
setwd("C:/Users/drmil/Dropbox/White House Lobbying/3YP/APSR Final/Replication Archive/visitor_logs_analyses")

load("tableSI11_cols1and2.RData")
load("tableSI11_cols3and4.RData")

################################################################################

# TO DETERMINE WHETHER THE DISTRIBUTIONS OF BETA COEFFICIENTS FROM THE HIGH-
# AND LOW-QUALITY ACCESS CONSTITUENT MODELS ARE STATISTICALLY DISTINGUISHABLE,
# WE USE THE hypothesis FUNCTION IN THE brms PACKAGE, WHICH DIFFERENCES THE
# COEFFICIENTS FROM EACH SAMPLING ITERATION AND ASSESSES THE DISTRIBUTION
# OF DIFFERENCES

# CLINTON

# LOBBYING EXPENDITURES
clinton_lobbyexp <- hypothesis(tableSI11_cols1and2, 
                               "anyvisitshighleq0_loglobbyexp_1lP1 = anyvisitslowleq0_loglobbyexp_1lP1")
clinton_lobbyexp

# CAMPAIGN CONTRIBUTIONS (NOTE: WE ASSESS THE SUM OF THE TWO RELEVANT COEFFICIENTS)

clinton_contribs <- hypothesis(tableSI11_cols1and2, 
                               "anyvisitshighleq0_made_donations_l1tol4 + anyvisitshighleq0_made_donations_l1tol4:logamt_donations_1ltol4P1 = anyvisitslowleq0_made_donations_l1tol4 + anyvisitslowleq0_made_donations_l1tol4:logamt_donations_1ltol4P1")
clinton_contribs

# INDUSTRY PARTISAN ALIGNMENT

# D vs. I

clinton_DvsI <- hypothesis(tableSI11_cols1and2, "anyvisitshighleq0_industry_pidI = anyvisitslowleq0_industry_pidI")
clinton_DvsI

# D vs. R

clinton_DvsR <- hypothesis(tableSI11_cols1and2, "anyvisitshighleq0_industry_pidR = anyvisitslowleq0_industry_pidR")
clinton_DvsR

# OBAMA

# LOBBYING EXPENDITURES

obama_lobbyexp <- hypothesis(tableSI11_cols3and4, 
                             "anyvisitshighleq0_loglobbyexp_1lP1 = anyvisitslowleq0_loglobbyexp_1lP1")
obama_lobbyexp

# CAMPAIGN CONTRIBUTIONS (NOTE: WE ASSESS THE SUM OF THE TWO RELEVANT COEFFICIENTS)

obama_contribs <- hypothesis(tableSI11_cols3and4, 
                               "anyvisitshighleq0_made_donations_l1tol8 + anyvisitshighleq0_made_donations_l1tol8:logamt_donations_1ltol8P1 = anyvisitslowleq0_made_donations_l1tol8 + anyvisitslowleq0_made_donations_l1tol8:logamt_donations_1ltol8P1")
obama_contribs

# INDUSTRY PARTISAN ALIGNMENT

# D vs. I

obama_DvsI <- hypothesis(tableSI11_cols3and4, "anyvisitshighleq0_industry_pidI = anyvisitslowleq0_industry_pidI")
obama_DvsI

# D vs. R

obama_DvsR <- hypothesis(tableSI11_cols3and4, "anyvisitshighleq0_industry_pidR = anyvisitslowleq0_industry_pidR")
obama_DvsR

################################################################################

# EXTRACT POINT ESTIMATES AND CI CUTOFFS FOR PLOTTING

beta_diffs <- data.frame("president"=c(rep("Clinton", 4),
                                       rep("Obama",4)), 
                         "covariate"=rep(c("LobbyExp", "Contribs", "D_vs_I", 
                                           "D_vs_R"),2), 
                         "pe"=c(clinton_lobbyexp$hypothesis$Estimate,
                                clinton_contribs$hypothesis$Estimate,
                                clinton_DvsI$hypothesis$Estimate,
                                clinton_DvsR$hypothesis$Estimate,
                                obama_lobbyexp$hypothesis$Estimate,
                                obama_contribs$hypothesis$Estimate,
                                obama_DvsI$hypothesis$Estimate,
                                obama_DvsR$hypothesis$Estimate), 
                         "ci_25"=c(quantile(clinton_lobbyexp$samples$H1, probs = c(0.025)),
                                   quantile(clinton_contribs$samples$H1, probs = c(0.025)),
                                   quantile(clinton_DvsI$samples$H1, probs = c(0.025)),
                                   quantile(clinton_DvsR$samples$H1, probs = c(0.025)),
                                   quantile(obama_lobbyexp$samples$H1, probs = c(0.025)),
                                   quantile(obama_contribs$samples$H1, probs = c(0.025)),
                                   quantile(obama_DvsI$samples$H1, probs = c(0.025)),
                                   quantile(obama_DvsR$samples$H1, probs = c(0.025))),
                         "ci_975"=c(quantile(clinton_lobbyexp$samples$H1, probs = c(0.975)),
                                    quantile(clinton_contribs$samples$H1, probs = c(0.975)),
                                    quantile(clinton_DvsI$samples$H1, probs = c(0.975)),
                                    quantile(clinton_DvsR$samples$H1, probs = c(0.975)),
                                    quantile(obama_lobbyexp$samples$H1, probs = c(0.975)),
                                    quantile(obama_contribs$samples$H1, probs = c(0.975)),
                                    quantile(obama_DvsI$samples$H1, probs = c(0.975)),
                                    quantile(obama_DvsR$samples$H1, probs = c(0.975))))

# RECREATING FIGURE 5

pdf(file = "fig5.pdf", family = "Times", height = 10, width = 12)
layout(matrix(c(1,2), ncol=1, byrow=TRUE),
       heights = c(0.85,0.15))
par(mar=c(5.1,18,2,.7))
plot(NULL, ylim=c(0.5, 16.5), xlim=c(-1.05,0.5), axes=FALSE,  bg = "black",
     tck=-.02, cex.axis=0.9, cex=1.5,
     xlab="", ylab="", yaxt="n", xaxt="n", main="")
abline(v=0, lty=2, col="gray")
points(beta_diffs$pe[1:4],
       c(15, 10, 5, 2),
       pch=19, cex=1)
points(beta_diffs$pe[5:8],
       c(14, 9, 4, 1),
       pch=17, cex=1)
segments(x0=beta_diffs$ci_25,
         x1=beta_diffs$ci_975,
         y0=c(c(15, 10, 5, 2), c(14, 9, 4, 1)), cex=1.5)
axis(1,at=c(-1.00, -0.50, 0.00, 0.50),
     labels=sprintf("%.2f",c(-1.00, -0.50, 0.00, 0.50)),cex.axis = 1.75, lwd=2)
axis(2,at=c(14.5, 9.5, 5, 4, 2, 1),
     las=2,
     labels=c(expression(~bold(~underline("Lobbying Expenditures"))), 
              expression(~bold(~underline("Campaign Contributions"))), 
              expression(~bold(~underline("Partisan Alignment"))),
              expression(~bold(~underline("(Republican)"))),
              expression(~bold(~underline("Partisan Alignment"))),
              expression(~bold(~underline("(Independent)")))
     ),
     tck=0,
     lwd = 0,
     line = 0,
     cex.axis=1.75)
mtext(TeX(r'($\beta_{High Quality}$ - $\beta_{Low Quality}$)'), side=1, line = 4, cex =2.25)
text(beta_diffs$pe[5:8],
     c(14, 9, 4, 1)+0.45,
     sprintf("%.2f",round(beta_diffs$pe[5:8], 2)), cex=1.5)
text(beta_diffs$pe[1:4],
     c(15, 10, 5, 2)+0.45,
     sprintf("%.2f",round(beta_diffs$pe[1:4], 2)), cex=1.5)
def_par <- par(no.readonly = T)
opar <- par(mar=c(0,0,0,0))
plot(0,0, type="n", axes=FALSE, xlab="", ylab="")
legend(-0.14,0.50,legend=c("Clinton",
                           "Obama"),
       pch=c(19,17), col=c("black"), ncol = 2,
       title="Presidential Administration", cex=1.75)
par(def_par)
dev.off()
